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Abstract 

Negative and positive transcriptional feedback loops are present in natural and synthetic genetic oscil- 
lators. A single gene with negative transcriptional feedback needs a time delay and sufficiently strong 
nonlinearity in the transmission of the feedback signal in order to produce biochemical rhythms. A single 
gene with only positive transcriptional feedback does not produce oscillations. Here, we demonstrate 
that this single-gene network in conjunction with a simple negative interaction can also easily produce 
rhythms. We examine a model comprised of two well-differentiated parts. The first is a positive feedback 
created by a protein that binds to the promoter of its own gene and activates the transcription. The 
second is a negative interaction in which a repressor molecule prevents this protein from binding to its 
promoter. A stochastic study shows that the system is robust to noise. A deterministic study identifies 
that the dynamics of the oscillator are mainly driven by two types of biomolecules: the protein, and the 
complex formed by the repressor and this protein. The main conclusion of this paper is that a simple 
and usual negative interaction, such as degradation, sequestration or inhibition, acting on the positive 
transcriptional feedback of a single gene is a sufficient condition to produce reliable oscillations. One gene 
is enough and the positive transcriptional feedback signal does not need to activate a second repressor 
gene. This means that at the genetic level an explicit negative feedback loop is not necessary. The 
model needs neither cooperative binding reactions nor the formation of protein multimers. Therefore, 
our findings could help to clarify the design principles of cellular clocks and constitute a new efficient tool 
for engineering synthetic genetic oscillators. 

Keywords: positive feedback - genetic oscillator - circadian clock - relaxation oscillator - hysteresis 

Abbreviations: NTF, negative transcriptional feedback; PTF, positive transcriptional feedback; QSSA, quasi- 
steady-state assumption 



Introduction 

Cellular clocks control important functions of the cell, such as circadian (24-hour) rhythms, cell cycle, 
metabolism and signaling. Clock operation appears to involve the coupling of two different types of oscilla- 
tors. The first are oscillators based on cytoplasmic reactions, such as phosphorylation [1] and oxidation [2,3]. 
The second are genetic oscillators depending on gene expression regulation [4,5]. In the last decade several 
synthetic genetic oscillators have been implemented in the laboratory [6-12]. The first mathematical model 
of a genetic oscillator was developed by Goodwin for periodic enzyme production [13]. This model was the 
groundwork for subsequent theoretical research on genetic oscillators in living systems, such as fungi and 
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Figure 1: Diagram of one-gene oscillators with negative and positive transcriptional feedbacks. A. 

Negative transcriptional feedback (NTF) created by a protein that represses the expression of its own gene. This 
NTF needs time delay and sufficiently strong nonlinearity in the feedback signal transmission in order to produce 
reliable oscillations. The time delay is created by intermediate reactions, such as the transcription and translation, 
reversibly phosphorylations or proteins shuttling between the nucleus and the cytoplasm. The nonlinearity can be 
created by reactions, such as protein cooperativity in the gene repression or formation of protein multimers. B. 
Positive transcriptional feedback (PTF) created by a protein that activates the expression of its own gene. This PTF 
needs a negative interaction in the feedback signal transmission in order to produce reliables oscillations. The negative 
interaction can be a degradation, sequestration, or inhibition carried out by a repressor molecule. 



flies [14-19]. In these models, the rhythms are generated by a gene with a negative transcriptional feedback 
(NTF) (Fig. 1A). This NTF needs time delay and sufficiently strong nonlinearity in the transmission of the 
feedback signal for preventing the steady-state stabilization of the system [20,21]. It has also been analyzed 
variants, involving two genes, of the model presented in the Fig. 1A [22]. 

Positive transcriptional feedbacks (PTFs) are also present in many cellular clocks [23-25]. Models with 
two or more genes involving PTFs have been studied in genetic oscillators [26-34] . In these models the PTFs 
increase the expression of repressor genes. It has been shown how PTFs produce bistability [35,36], increase 
the robustness of cellular clocks [37, 38] and could provide robust adaptation to environmental cycles [39] . 
Previously, it has been demonstrated that a single gene with only PTF does not produce oscillations [40]. 
Here we study a model with a simple condition to produce biochemical rhythms in a single gene with PTF 
(Fig. IB). We chose a circadian period for the oscillator due to its relevance in biological systems. This model 
is based on two common features of genetic oscillators [4,21,26,28,38]. The first is a PTF created by a protein 
that activates the transcription of its own gene. The second is a negative interaction in which a repressor 
inhibits the activity of this protein. We performed stochastic and deterministic simulations that yielded 
similar results. The stochastic simulations show that the genetic oscillator is robust to noise. This noise is 
introduced in living cells by the stochasticity of gene expression [41,42]. By means of a reduced deterministic 
model, we show that the oscillations exhibit limit-cycle behavior. This means that if a disturbance is applied 
to the system, the oscillations return to the original periodic solution [43,44]. Also we show that this biological 
clock can be classified as a relaxation oscillator [28,43,44]. This type of clock is sometimes called hysteresis 
oscillator [26,45] or amplified negative feedback oscillator [21,25]. The relaxation oscillator comprises fast 
and slow oscillation creation stages. In our model these oscillations are characterized by sawtooth waveforms. 
Finally, we explain how the negative interaction works through a comparison with the dynamics of the typical 
enzymatic reaction. We show that the rate of the negative interaction is amplified by the PTF and has a 
saturation point. 
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Figure 2: Model of a one-gene oscillator with PTF. The model is composed of two well-differentiated parts. 
The first part is a positive feedback loop in which a gene (G) is transcribed into mRNA (M). In turn, M is translated 
into protein (A). This protein is a transcription factor of its own gene and increases the transcription rate when it 
binds to the promoter. The positive feedback needs a second part, consisting of a negative interaction in order to 
obtain reliable oscillations. In this part repressor molecules (R) enter the system at a constant rate. R inhibits the 
function of A. Specifically, R binds to A and forms the complex C. In this complex, A is not able to bind to its 
promoter. R is not degraded together with A and can be used several times. Therefore, R can be thought of as a 
protease, a protein that sequesters A or any other molecule that binds to and inhibits the function of A as explained 
above. The zigzag arrows stand for degradations. A different version of the model can be formulated with the negative 
interaction acting over M instead of over A. 

Results 

Model and simulations 

The model is a simple one-gene network with two well-differentiated parts (Fig. 2). The first is a PTF created 
by a protein A, which is a transcription factor of its own gene. When this protein binds to its promoter 
the transcription rate increases. The second part is a negative interaction in which a repressor molecule R 
prevents A from binding to its promoter. The molecule R can be thought of as a protease, as a protein that 
sequesters A, or as any other molecule that inhibits the function of A as shown in Fig. 2. A different version 
of the model can be formulated in which the negative interaction acts on the mRNA molecules instead of on 
protein A. 

Eleven biochemical reactions provide a full description of the model (see (3) in the section Methods: 
Biochemical reactions and rates). The system is assumed to have a uniform mixture of biomolecules. For 
this reason, we did not take into account diffusion processes. In this approach, the dynamics of the biochemical 
reactions (3) can be described by two different formalisms known as stochastic and deterministic approaches 
(see Methods: Deteministic and stochastic simulations for more details). These two approaches can lead 
to different behaviors. The stochastic dynamics of the reactions (3) were simulated using the Gillespie 
algorithm [46] and the deterministic dynamics using the following ordinary differential equations: 



where the variables and rates are described in the section Methods: Biochemical reactions and rates. We 
used standard values within the diffusion limit for the rates [18,38,47]. 
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Figure 3: Stochastic and deterministic simulations of the model. A, C, E, G, I, K. Stochastic time 
evolution of the protein (A), repressor (R), protein-repressor complex (C), mRNA (M), gene (G) and activated gene 
(G a ), respectively. B, D, F, H, J, L. Deterministic time evolution of A, R, C, M, G, and G a , respectively. In 
both simulations, the time evolution of A (A and B), R (C and D) and C (E and F) are very similar except for 
the presence of fluctuations in the stochastic case. This phenomenon is more pronounced in the time evolution of M 
(G and H). The oscillations in C show sawtooth waveforms. There is a single gene in the model; hence G + G a = 1 
molecule. In the time evolution of G (I and J) and G a (K and L), the stochastic simulation shows discrete transitions 
between and 1 molecules. By contrast, the deterministic simulation shows unrealistic continuous transitions. The 
time evolution of G a shows that the gene is activated most of the time (K and L). 



The stochastic approach is more realistic than the deterministic simulation because it takes into account 
the randomness of the chemical reactions. This randomness produces fluctuations in the number of molecules. 
We fitted the reaction rates to obtain circadian oscillations in the stochastic simulation. Then, we compared 
the results with the deterministic simulation (Fig. 3). For both simulations the time evolution of the protein 
(A), repressor (R), protein-repressor complex (C) and mRNA (M) arc very similar. The main difference 
is the appearance of fluctuations in the stochastic case around the number of molecules predicted by the 
deterministic approach. The fluctuations are more evident in the time evolution of M (Fig. 3G) than in 
the other biomolccules. This is because the number of M molecules oscillates in a lower range than A, 
R and C. The oscillations in C are characterized by sawtooth waveforms. On the other hand, there are 
differences between the stochastic and deterministic time evolution of the gene. There is a single gene in the 
model, which can be deactivated (G) or activated (G Q ). Therefore, G + G a = 1 molecule. The stochastic 
simulation shows realistic discrete transitions between and 1 molecules (Figs. 31 and 3K). By contrast the 
deterministic simulation shows unrealistic continuum transitions (Figs. 3 J and 3L). In both cases, however, 
the qualitative behavior is the same. Most of the time the gene is activated by A, although it is deactivated 
for a short time when the number of A in the oscillations is low. 

Model robustness to noise 

The fluctuations in the stochastic simulation are the source of so-called intrinsic noise [41,42]. In the genetic 
oscillator, this intrinsic noise generates variability in both the amplitude and period of the oscillations. The 
phase plane defined by C and A illustrates this variability very clearly (Fig. 4A). The deterministic phase 
plane is a well-defined curve because the oscillations are identical (dashed line in Fig. 4A). In contrast, the 
stochastic phase plane is a curve that spreads around the deterministic curve due to intrinsic noise (solid line 
in Fig. 4A). We used the amplitude and period histograms, and the autocorrelation function to quantify the 
effect of this intrinsic noise on A oscillations. The results are similar to circadian models with more chemical 
reactions [18,26]. The amplitude histogram shows a mean of 6,723 molecules and a standard deviation of 
858 molecules (Fig. 4B). The period histogram shows a mean of 24.3 hours and a standard deviation of 1.7 
hours (Fig. 4C). In contrast, the absence of intrinsic noise in the deterministic simulation produces identical 
A oscillations with lower amplitude and period equal to 6,164 molecules and 23.6 hours, respectively. On the 
other hand, the autocorrelation function shows a half-life time of about 120 hours (Fig. 4D). 
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Figure 4: Model robustness to noise. A. Stochastic and deterministic phase plane defined by C and A. The 

general shape of the curves is the same in both cases, but the deterministic curve (dashed line) is well defined because 
the oscillations are identical. In contrast, the stochastic curve (solid line) spreads around the deterministic one. B, C. 
Amplitude and period histograms of the stochastic simulation of A, respectively. D. Autocorrelation of the stochastic 
oscillations in the number of A molecules. The half-life of the autocorrelation is about 120 hours (intersection of 
dashed lines). E-H. Model robustness to intrinsic noise when the number of molecules is low. The changed rates are 
k± = 1000 hour -1 , fcs = 5000 hour" 1 , k7 = 25.5 molecules -1 hour -1 , k% = 132.6 hour -1 and kg = 1 molecules hour -1 . 
In particular, we multiplied the rates ki and fcs by 100 to obtain a low number of M molecules. Simultaneously, we 
multiplied the rates k*i and fcg, and divided the rate fcg by 51 to obtain a low number of R and C molecules. The 
initial conditions are G a o = 1 and Go = Mq — Aa = Ro — Co = molecules. The mean value of M is 0.48 molecules. 
I. Model robustness to extrinsic noise. Scatter plot of amplitude versus period that shows the robustness of the model 
to parameter variation (data is presented in Table SI). Two stochastic simulations were performed for each parameter 
in which the value was increased and decreased by 15%. The x and y coordinates of each data point correspond to 
the mean values of the period and amplitude, respectively. The horizontal and vertical error bars are the standard 
deviation of the period and amplitude, respectively. The intersection between dashed lines shows the point obtained 
without changing the value of any rate (Figs. 4B and 4C). (B, C, D and each data point in I were calculated for 
1,000 successive cycles. We assumed that a cycle occurs if the number of proteins A increases to 1,000 molecules and 
then decreases to 700 molecules. The amplitude was calculated as the greatest number of A molecules in each cycle. 
The period was calculated as the time interval that it takes the number of proteins A to reach 1,000 molecules for the 
first time in two successive cycles.) 
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The stochastic approach produces good oscillations in A even when there are fewer than 30 molecules of 
M, R and C. (Figs. 4E-H). We changed the value of some rates to obtain this simulation as in ref. 38 (see 
caption of Fig. 4). In the deterministic approach, where intrinsic noise is not present, these changes do not 
alter the dynamics of A significantly and produce a low number of M , R, and C molecules. In particular, the 
amplitude and the period are slightly lower (Fig. SI). In the stochastic simulation the rate changes reduce 
the amplitude and period means to 6,166 molecules and 21.3 hours, respectively (Fig. S2). The effects of 
intrinsic noise is now more pronounced because the number of M, R, and C molecules is low. This is reflected 
in an increase of the amplitude and period standard deviations to 2,132 molecules and 5.2 hours, respectively 
(Fig. S2). 

In cells, there are also fluctuations in the number (or activity) of molecules such as polymerases, ribosomes 
and degradation machinery. These fluctuations are the source of so-called extrinsic noise [41,42]. We 
performed stochastic simulations varying the parameters in order to account for some aspect of extrinsic 
noise in the robustness study of the model. The results show that this oscillator is robust to small parameter 
variations (Fig. 41) like more other complex models of genetic oscillators [27]. The largest amplitude and 
period changes occurred for variations in k% (sec Table SI). The changes in the mean period and amplitude 
were always less than 15% and 31%, respectively. Particularly, variations in the rates fci, k—i, k 2l k§, k*? and 
fcio produced changes of less than 3% and 8% in the mean period and amplitude, respectively. The changes 
in the standard deviation of the period and the amplitude were always less than 13% and 27%, respectively. 



Reduced deterministic model 

To identify the types of biomoleculcs mainly responsible for oscillations, it is useful to reduce the deterministic 
model by means of the quasi-steady-state assumption (QSSA) [43,48]. This approximation differentiates 
between fast and slow variables. The greater the time-scale separation between the variables the more 
accurate the approximation is. In this approach it is assumed that fast variables quickly reach the equilibrium, 
i.e., their derivatives are zero. This assumption means that slow variables are responsible for the system 
dynamics. In this model, we assumed that the fast variables are G, G a , M and R, and the slow variables are 
A and C. Then, the set of Eq. (1) can be simplified to 



dA/dt = i±ia >+y - ^ 

7+A S+A 



dC/dt 



k 9 A - Sk 8 C 

sTa ' 



(2) 



where a = Gtk-ik 2 k 5 /kik4, (3 = Gtk^k^/ki, 7 = k-i/ki, 5 = ki /kr and G t = G + G a . A good way to check 
if this approximation is correct is to compare the numerical solution of the complete and the reduced systems. 
Both numerical solutions agree except for quantitative differences in the period and the amplitude (Figs. 5A 
and 5B). These differences arc due to the fact that the time-scale separation between fast and slow variables 
is not large enough for QSSA to be more accurate. Despite these differences, we can conclude that A and C 
are mainly responsible for the system dynamics. The other types of biomoleculcs can be considered to be at 
equilibrium. The fluctuations in the fast variables do not significantly affect the system dynamics [38]. This 
explains the robustness of the model when the number of molecules is low (Figs. 4E-H). In fact, the system 
produces reliable oscillations even if the average of M is less than one molecule (Fig. 4H), and, surprisingly, 
even when the driven C molecules oscillate in a range of less than 30 molecules (Fig. 4F) . 

The oscillations in the reduced deterministic model exhibit limit-cycle behavior (thin solid line in Fig. 
5C). Therefore, if an external disturbance is applied to the oscillator, the system will go back to oscillating 
with the period and amplitude of its limit cycle. The unstable fixed point of the system is Co = 552.4 and 
Aq = 56.3 molecules (circle in Fig. 5C). For a bifurcation analysis of parameters ks and fcg indicating the 
range of values that produces limit-cycle oscillations, see Methods: Bifurcation diagram. 

This genetic clock belongs to the so-called relaxation oscillators [28,43,44]. The mechanism responsible 
for the oscillations is represented by the nullclincs An and Cn (Fig. 5C). These nullclines are the solution 
of the equations dA/dt = and dC/dt = 0, respectively. The nullcline Cn is a straight line and the nullclinc 
An has the characteristic "Z" shape of relaxation oscillators [43-45]. The shape of the A nullcline is the same 
as the hysteresis diagram obtained if C is assumed constant (Fig. S6). Therefore, this genetic clock contains 
some features of hysteresis in its oscillatory mechanism. The A nullcline has two branches that we can call 
"high" and "low" (Fig. 5C). These branches are steady states if the C is a constant (Fig. S6). In each 
oscillation the system switches from one branch to the other using the number of C molecules as a transient 
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Figure 5: Model of the genetic oscillator reduced by QSSA. A, B. Comparison between the reduced (solid 
line) and complete (dashed line) deterministic simulation of the time evolution of A and C, respectively. C. Phase 
plane. Limit cycle (thin solid line) and nullclines An (thick solid lines) and Cjv (thick dashed line). The unstable 
fixed point of the system (marked by circle o) is Co = 552.4 and Ao = 56.3 molecules. The nullclines An and Cjv are 
the solution of equations dA/dt = and dC/dt = 0, respectively. The two branches in the nullcline An are called 
"high" and "low". D. Slow and fast stages in the reduced system. The solid line is C and the dotted line is A. C 
exhibits a sawtooth waveform. (The arrows in C and D represent the direction of the oscillations. One and two 
arrows mean slow and fast stages, respectively). 

signal. This process can be explained following the limit-cycle trajectory. When A and C are about 1 and 
200 molecules, respectively, their number increases until A reaches its maximum of about 7,330 molecules 
and C reaches about 650 molecules. This is the transient from the low to the high branch. Then, the number 
of A molecules is reduced to about molecules, whereas C reaches its maximum of about 1,260 molecules. 
This is the transient from the high to the low branch. Finally, the number of C molecules is quickly reduced 
and the trajectory moves along the nullcline An, returning to the starting point where a new cycle begins. 

This genetic clock is characterized by containing fast and slow stages. The time evolution of C shows 
these two well-differentiated stages (Fig. 5D). In the slow stage A ^> S and kgA S> 6k$C, then the second 
differential equation in (2) can be approximated by dC/dt ks kg. In this stage, therefore, the number of C 
molecules increases linearly according to equation C oc kgt. In the fast stage A <C S and kgA <C <5fcgC, then 
the second differential equation in (2) can be approximated by dC/dt ~ — ksC. In this stage, the number of 
C molecules decays exponentially according to equation C cx cxp(— kgt). The two stages play different roles. 
The slow stage is characterized by the formation of a pulse of A molecules. On the other hand, the decay 
of C into R in the fast stage provides the necessary conditions for a new pulse. These two stages produce 
oscillations in C with sawtooth waveforms (solid line Fig. 5D). 
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How the negative interaction works 



The negative interaction decreases the number of free A molecules and takes the system back to the start 
of a new cycle. The detailed explanation of how this interaction works is related to the dynamics of the 
typical enzymatic reaction S + E N D — » P + E, where S, E, D and P are the substrate, enzyme, com- 

C_l 

plex substrate-enzyme and product, respectively. The total number of enzymes (E t = E + D) is constant 
in the system. The rate of catalysis in this reaction is defined as v = dP/dt = C2D. The value of this rate 
can be approximated by QSSA. The result of this approximation is the well-known Michaelis-Menten equa- 
tion v « V max S/(KM + S), where V max = c 2 E t and Km = (c_i + c 2 )/ci [43]. In this equation, the rate v 
increases asymptotically as a function of 5". The rate v reaches a maximum value {V max ) when the amount 
of S is large compared with the constant Km- In this situation, the enzymes are saturated because most are 
part of complex D, and adding more S does not increase the rate v. Therefore, D « E t , and the rate of the 
catalysis v reaches the constant value C2E t . 
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Figure 6: Rate of the negative interaction. A. Rate of the negative interaction (v = kgC). This rate represents 
the number of degraded A molecules per hour. The graph was plotted by multiplying the number of C molecules 
in Fig. 3F by k$. The circle (o) indicates the saturation point. At the saturation point the rate increases linearly 
(v oc kskgt) because new R molecules enter the system at rate kg. The square (□) indicates the maximum rate of 
the negative interaction {v max — 3,180 molecules/hour). B. Plot of v max against ku, where ku is the rate of the 
reaction: C — > (f>. Each point corresponds to a deterministic simulation with ku equal to 0.0, 0.1, 0.2, 0.3, 0.4, 0.5 
and 0.6 hour -1 , respectively (see Fig. S5 for more detailed information). The oscillations stop when fen = 0.6 hour -1 
(o) (Fig. S5G). 

In our model, the negative interaction is A + R s C — > R, where we assumed fc_7 = to simplify 

fc— 7 

the model. We can think of A, R, and C as S, E and D, respectively. Therefore, the rate of the negative 
interaction can be defined as v = fcgC (Fig. 6A). This rate represents the number of degraded A molecules 
per hour. The negative interaction works as follows. The number of A molecules increases quickly due to the 
positive feedback. This rise causes most of the R molecules to bind to A molecules forming the complex C. At 
this point, the system reaches the saturation level (circle in Fig. 6A). The total number of repressor molecules 
in the system is Rt = R + C. Therefore, at the saturation point, C sa Rt and the rate v reaches the value 
kgR t . The negative interaction is not fast enough to decrease the growth of A molecules immediately after 
the saturation point is reached. This is because the number of Rt molecules is low at this point. Nevertheless, 
new R molecules enter the system at rate kg. Therefore, Rt increases linearly over time (Rt cx kgt) compared 
with the enzymatic reaction in which E t is constant. This means that the rate of the negative interaction 
increases linearly according to equation v oc kgkgt. The value of v increases until the negative interaction is 
fast enough to reduce the number of A molecules and take the system back to the start of a new cycle. The 
maximum rate reached by the negative interaction is v max = 3,180 molecules/hour (square in Fig. 6A). 

In this model there is not an explicit negative feedback loop at the genetic level. It has been conjectured 
that all biochemical oscillators involve some sort of negative feedback loop [21]. In this genetic clock, an 
effective negative feedback loop appears in the reduced model (see the section Methods: The Jacobian matrix). 
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Intuitively, this effective negative feedback loop can be explained as follows: when C is rare, A is increased by 
the positive feedback. This rise in the production of A leads to the accumulation of C, which in turn increases 
v. This accumulation of C increases until the negative interaction is fast enough to reduce the number of A 
molecules. In this model, we assumed that C is not degraded. If this complex is degraded according to the 

reaction C — ^ 4>, v increases at a slower rate, and its maximum value (v max ) is reduced (Figs. 6B and S5). 
The oscillations stop when fen = 0.6 hour -1 (Fig. S5G), because not enough C is accumulated in order to 
increase v. 

This genetic oscillator does not need cooperative binding reactions nor the formation of protein multimcrs, 
in contrast to the one-gene oscillator with TNF (Fig. 1A). It has been demonstrated that protein sequestration 
produces an effective high nonlincarity [49,50]. But this high nonlinearity is not observed if the repressor 
molecule is recycled [49]. In our model the repressor R can be used several times. Therefore, the negative 
interaction does not produce an effective high nonlinearity (see Supporting Information: Text SI). 

Discussion 

Genetic networks with NTFs and PTFs play an important role in cellular clocks. In this paper, we provided 
a simple model illustrating that a single gene with PTF has also the potential to produce reliable oscillations. 
The sufficient additional requirement is a simple and usual negative interaction of degradation, sequestration 
or inhibition acting on the positive feedback signal. The model presented in this article has a different 
oscillatory mechanism than the well-established NTF one-gene oscillator model. Our model can be classified 
as a relaxation oscillator. A two-gene model has been proposed as a different way of producing reliable 
circadian oscillations in cellular clocks [26], which also is a relaxation oscillator. This two-gene model is 
important because it is robust to noise [38] . The model introduced in this paper is a simpler way to produce 
relaxation oscillations than the previous two-gene oscillator. A comparison with our model reveals that the 
activation of the repressor gene is not a necessary condition to produce reliable circadian oscillations in the 
two-gene oscillator. We demonstrated that our model produces circadian oscillations that are just as robust to 
noise as the two-gene oscillator and other more complex models [18,27]. Similarly to the two-gene oscillator, 
our model produces good oscillations when the average number of mRNA molecules is less than one. In fact, 
the number of proteins oscillates satisfactorily even when the other types of molecules involved in the clock 
are less than 30. Therefore, this model is a simpler genetic relaxation oscillator than the current two-gene 
clocks [25]. Our model does not need the activation of a second repressor gene by the PTF, cooperative 
binding reactions nor the formation of protein multimers. 

A single gene with PTF and a negative interaction in the feedback signal is an alternative and simple way 
of generating reliable oscillations. Our study suggests that PTF, besides increasing robustness in cellular 
clocks, could be more directly and deeply involved in the production of oscillations than at first thought. 
Further research is necessary to elucidate the presence and the role of this genetic oscillator in natural 
cellular clocks. On the other hand, thanks to its simplicity, this model has the potential to be a new tool for 
engineering synthetic genetic oscillators. In this case the period and amplitude of the oscillations could be 
possibly controlled by externally manipulating the entry rate of the repressor molecules. 
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Methods 



Biochemical reactions and rates 

The biochemical reactions that fully describe the model in the Fig. 2 are as follows: 
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Complex creation: 


R 




C 


Complex decay into R: 


C 


^R 




R creation (or entry): 


<t>- 


k ^R 




R degradation (or exit): 


R 


— > 9, 





(3) 



where G denotes the gene without A bound to its promoter, M denotes mRNA transcribed from G, A 
denotes the activator protein translated from M, G a denotes the gene with A bound to its promoter, R 
denotes the repressor and C denotes R bound to A. All the biochemical species are measured in molecules. 
The description of the rates is as follows: k\ is the binding rate of A to the promoter of G, is the unbinding 
rate of A from the promoter of G, ki is the basal transcription rate, k^ is the activated transcription rate, k^ 
is the degradation rate of M, k$ is the translation rate, k 6 is the degradation rate of A, kj is the binding rate 
of R to A, k$ is the decay rate of C into R, kg is the creation (or entry) rate of R and fcio is the degradation 
(or exit) rate of R. 

We used standard values within the diffusion limit for the rates [18,38,47]. They are as follows: ki = 1 
molecules -1 hour -1 , k-i = 50 hour -1 , ki = 50 hour -1 , k% = 500 hour -1 , k& = 10 hour -1 , k§ = 50 hour -1 , 
fcg = 0.1 hour -1 , k 7 = 0.5 molecules -1 hour -1 , fc 8 = 2.6 hour -1 , fc 9 = 51 molecules hour -1 and fcio = 1 
hour -1 . The cell has a single copy of the gene: Gt = G + G a = 1 molecule. The initial conditions are: 
Go = 0, G a o = 1, Mq = 5, Aq = 1000, Ro = 5, and Go = 1200 molecules. The initial conditions have been 
chosen to obtain a first cycle with an amplitude similar to the limit-cycle oscillations. Note that the rates k\ 
and kj include the volume of the system V. Hence, these rates can be written as k\ = k\jV and kj = k?/V, 
where the rates fcj and kj are expressed in M hour -1 . In order to generate circadian oscillations, first, we 
varied all the reaction rates, according to the values used in refs. 18, 38 and 47, until we got oscillations with 
a period of around 24 hours in the stochastic simulation. Then we fine-tuned the oscillations varying rates 
/c8 and fcg until a period closer to 24 hours was achieved. 

Deteministic and stochastic simulations 

Models based on chemical reactions in a well stirred system are usually described by two different formalisms 
from a mathematical point of view: 

Deterministic: this formalism is suitable for large numbers of molecules. It is described by a set of coupled 
ordinary differential equations that follow the law of mass action. These equations are called reaction rate 
equations and they can only be solved analytically for simple systems. For more complex systems numerical 
methods are necessary. In this approach the amount of each chemical species and the time are continuous. 
The velocity at which reactions occur is given by the reaction rate constants k, or simply rate. 

Stochastic: this formalism is suitable for small numbers of molecules because it takes into account the 
randomness of the chemical reactions. It is described by the so-called master equation, which is the time 
evolution of the probability that the system has a certain number of molecules of each chemical species at 
time t. Few systems can be solved analytically with the master equation. It is possible, however, to simulate 
the stochastic behaviour with the Gillespie algorithm [46]. In this approach the amount of each chemical 
species and the time are discrete, and the rates k turn into probabilities. 
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Bifurcation diagram 



We calculated the bifurcation diagram for parameters kg and kg. These are key parameters for two reasons. 
First, the rate of the negative interaction v is proportional to kg and kg when the saturation point is reached. 
Second, the fast and slow stages in the relaxation oscillations depend on fcg and kg, respectively. Specifically, 
wc studied the range values of kg that produce stable oscillations through a bifurcation diagram. Then we 
studied how this range changes when the parameter kg varies. 

The bifurcation diagram of the reduced model depending on kg shows two Hopf bifurcation points (Fig. 
S3A). The first Hopf bifurcation appears at kg — 4.78 molecules hour -1 and the second at kg = 217.6 
molecules hour -1 . Most of the values of kg between these two points produce stable oscillations. Only 
for a short range of values around these points are the oscillations unstable (white circles in Fig. S3 A). 
The oscillations have an amplitude of from 2,000 to 16,000 molecules, and a period of from 7 to 170 hours 

kg 

(Fig. S3B). The velocity of the reaction <f> — > R in (3) does not depend on any biomolecule involved in the 
oscillator. Therefore, parameter kg can be interpreted as an external signal controlling the behaviour of the 
clock. 

The variation of parameter kg changes the position of the two Hopf bifurcation points (white circles in 
Fig. S4). The different positions of these points define the regions with stable oscillations depending on the 
values of kg, and kg (regions I and II in Fig. S4). If parameter kg, is increased, the range of values of kg that 
produces stable oscillations decreases. This range shrinks faster if kg is greater than 20 hour -1 . We plotted 
an equivalent graph for the stochastic model because it is more realistic than the reduced graph (black circles 
in Fig. S4). In particular, we assumed that oscillations occurs in a region if the correlation in the first period 
is greater than 0.2. The stochastic model produces oscillations in the regions II and III (Fig. S4). The range 
of oscillations in the complete deterministic model is close to the region II. 



The Jacobian matrix 



The Jacobian matrix of the reduced system (2) is: 



J = 



an a 12 

0-21 a 22 



( 7j3-a S(k 9 + kgC) 

(7 + A) 2 (S + A) 2 
S(k 9 + kgC) 



\ 



{5 + Af 



kgA \ 

'5 + A 
8kg 

'JTaJ 



(4) 



where the clement a\2 and 022 are always negative, the element 021 is always positive and the element an 
can be positive or negative depending on the values of the rates. With the rates given in the section Methods: 
Biochemical reactions and rates and the fixed point of the reduced system (Fig. 5C) the sign pattern for the 
Jacobian matrix is: 



J 



(5) 



A two-component negative feedback loop is created in the reduced model because 012021 < (sec Chapter 9 
of the reference [48]). The Jacobian matrix (5) has a tipically sign pattern that produces Hopf bifurcation 
in chemical systems with two variables [43,48]. The two-component systems with this sign pattern in the 
Jacobian matrix arc called activator- inhibitor models [481. 



Software 

Code for stochastic and deterministic simulations was written in FORTRAN and XPPAUT (http://www. 
math.pitt . edu/~bard/xpp/xpp .html), respectively. Simulations have been contrasted using CAIN software 
(http://cain.sourceforge.net/). The stability analysis to determine steady states and limit cycles was 
performed with XPPAUT. The histograms and autocorrelation function were plotted using FORTRAN and 
GNU Octave (http://www.gnu.org/software/octave/). The code for complete and reduced deterministic 
simulations in XPPAUT is available in File SI and File S2. The code for stochastic and deterministic 
simulations in CAIN is available in File S3. 
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Supporting Information: Figures 
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Figure SI. Time evolution of A with and without a low number molecules. Comparison between 
deterministic simulation of the time evolution of A with (dashed line) and without (solid line) a low number 
of M, R, and C molecules. (Solid line graph: the values of the parameters are as in the section Methods: 
Biochemical reactions and rates. Dashed line graph: the changed rates are k^ = 1000 hour , fcs = 5000 
hour -1 , k 7 = 25.5 molecules -1 hour -1 , fc 8 = 132.6 hour -1 and kg = 1 molecules hour -1 .) 
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Figure S2. Amplitude and period histograms of the stochastic simulation of A. A, B. Amplitude 
and period histograms of the stochastic simulation of A, respectively. The values of the parameters are as in 
the section Methods: Biochemical reactions and rates but now we set = 1000 hour -1 , fcg = 5000 hour , 
kj = 25.5 molecules -1 hour -1 , k% = 132.6 hour -1 and kg = 1 molecules hour -1 . (A and B were calculated 
for 1,000 successive cycles. We assumed that a cycle occurs if the number of proteins A increases to 1,000 
molecules and then decreases to 700 molecules. The amplitude was calculated as the greatest number of A 
molecules in each cycle. The period was calculated as the time interval that it takes the numbers of proteins 
A to reach 1,000 molecules for the first time in two successive cycles.) 
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Figure S3. Bifurcation diagram of the reduced model. A. Bifurcation diagram depending on kg. The 
solid/dashed line represents stable/unstable fixed points. Black/white circles are the maximum and minimum 
values of A during unstable/stable oscillations. HB denotes a Hopf Bifurcation point. HBi and HB2 appear 
when the value of kg is 4.78 and 217.6 molecules hour -1 , respectively. B. Period of the stable oscillations in A. 
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Figure S4. Oscillatory regions in the reduced and stochastic models depending on fc§ and kg. 
Region I. Oscillations in reduced model. Region II. Oscillations in both reduced and stochastic model. 
Region III. Oscillations in the stochastic model. Region IV. No oscillations in any model. White circles 
represent the locus of Hopf bifurcations in the reduced model (data are presented in Table S2). Black circles 
represent locus of oscillations in the stochastic simulation (data are presented in Table S3). We assumed in 
the stochastic case that oscillations occur in a region if the correlation in the first period is greater than 0.2. 
(The lines connecting circles are designed to clearly single out the different regions.) 
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Figure S5. Rate of the negative interaction for different values of fen. Rate of the negative interac- 
tion (v = kgC) for different values of fen, where kn is the rate of reaction C — > <fi. Deterministic simulations 
A, B, C, D, E, F and G correspond to fen equals 0.0, 0.1, 0.2, 0.3, 0.4, 0.5 and 0.6 hour -1 , respectively. 
The values of the other parameters are as in the section Methods: Biochemical reactions and rates. The 
oscillations stop when kn = 0.6 hour -1 (G). If kn is increased, v increases slower, and its maximum value 
(v max ) is lower. The value of v max corresponds to the peak of the oscillations (v max is the value of the steady 
state in G). 



19 




j» 100 ; 




200 600 1000 1400 
C (molecules) 

Figure S6. Hysteresis diagram. Hysteresis diagram depending on C . The curve is the solution of the 
equation dA/dt = 0, where C is assumed constant. The two solid lines in the diagram are the two stable 
steady states "high" and "low" as a function of C. The dashed line represents the unstable points in the 
diagram. 



20 



Supporting Information: Tables 

Table SI. Data points of Fig. 41. 





Period of A 


Amplitude of A 


Changed rate 


Mean (hours) 


S.D. (hours) 


Mean (molecules) 


S.D. (molecules) 


none 


24.3 


1.7 


6723 
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1 i r 7 

l.lofei 


23.9 


1.8 


6554 
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0.85fci 


24.6 


1.6 


6792 
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i i r 7 

1.15fe_i 


24.5 


1.6 
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0.85/c_i 
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1.9 
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1.15K2 
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0.85fc 2 
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865 


i i r 7 
1.15fc 3 
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0.85fc 3 
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1 1 P* 7 

1.15K4 


21.9 


1.6 


5263 
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0. 85/1-4 


27.2 


1.8 
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i i r 7 

1.15k 5 
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Table S2. Data points of locus Hopf bifurcation in reduced model (Fig. S4). 



kg (hour 1 ) 


fcmzn ( mo l ecu l es hour X ) 


k™ ax (molecules hour 1 ) 


0.1 
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220.3 
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4.70 
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2.6 


4.78 


217.6 


3 


4.80 


217.1 


3.5 


4.82 
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5 


4.90 


214.9 


10 


5.17 


209.6 


20 


5.75 


199.3 


30 


6.39 


189.7 


40 


7.10 


180.8 


80 


10.5 
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120 


14.8 


127.1 


160 


20.1 


108.4 


200 


26.9 


91.9 
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81.5 
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41.4 


68.8 
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44.2 


65.3 
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48.3 


60.5 


261 


49.6 


59.1 


262 


51.5 


57.0 


262.6 


53.7 


54.8 
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Table S3. Data points of locus of oscillations with less than 20% of correlation in the first 
period in the stochastic model (Fig. S4). 



k$ (hour 1 ) 


k™ m (molecules hour x ) 


fornax ( mo l ecu l es hour 1 ) 
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16 
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Supporting Information: Text 

Text SI. The negative interaction does not produce an effective high nonlinearity. 

ft has been demonstrated that protein sequestration produces an effective high nonlinearity [49, 50]. But 
this high nonlinearity is not observed if the repressor molecule is recycled (see equation S9 and figure S5 in 
[49]). The biochemical reactions that describe the negative interaction are as follows: 



A creation (or entry): 




A degradation: 


.4 


<P 


Complex creation: 


R 


+ A^ 


Complex decay into R: 


C 


%R 


R creation (or entry): 


cf>- 


*%R 


R degradation (or exit): 


R 





c 



k 7 k 9 + k 6 k 10 

h 

_ fk 7 kg 



(k 7 k 9 + k 6 k w )k 8 ' 

where we observe no nonlinearity in output A as a function of input flux /. 



(6) 



The dynamics of these reactions are described by the following EDOs: 

dA/dt = f - k 6 A - k 7 AR 

dR/dt = -k 7 AR + k 8 C + k 9 - k 10 R (7) 
dC/dt = k 7 AR - k 8 C, 

As in [49], these equations can be solved at steady state to yield: 

fklQ 



A 

kq , , 

R=T L (8) 
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Supporting Information: Files 

File SI. Complete deterministic model (XPPAUT software). 

Available in: 

http : / /www .plosone . org/article/f etchSingleRepresentation .action?uri=inf o :doi/10 . 1371/ journal . 
pone. 0027414. sOll 

File S2. Reduced deterministic model (XPPAUT software). 

Available in: 

http : //www .plosone .org/article/f etchSingleRepresentation .act ion?uri=info :doi/10 . 1371/ journal . 
pone. 0027414. s012 

File S3. Stochastic and deterministic model (CAIN software). 

Available in: 

http : / /www .plosone .org/article/f etchSingleRepresentation .act ion?uri=inf o :doi/10 . 1371/ journal . 
pone. 0027414. s013 
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